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Abstract 


This paper develops a beam pattern design procedure for general multidimensional 
irregular sonar arrays that incorporates the not well understood effects of array geom- 
etry into the design process. The procedure is implemented by generating a “penalty 
function” in a spectral covariance function form. Processing the penalty function causes 
beam pattern high sidelobes to be penalized and the main lobe to be emphasized. This 
is accomplished by forming the penalty function in terms of an isotropic noise field of 
specified strength modified with a finite sector of low coherent energy and stabilized with 
incoherent sensor noise. By inputting the penalty function into a minimum variance 
beamformer, the beam pattern and aperture weights are calculated based on the given 
array geometry. The beamformer used is Capon’s Maximum Likelihood Method. The 
array used to test the procedure is located on a sixty degree sector of a cylindrical sur- 
face. The procedure is implemented by two different methods, each with some desirable 
characteristics. One method suppresses sidelobes directly by the placement of nulls. The 
other method suppresses sidelobes indirectly by the enhancement of the main lobe with 
anti-nulls. Both methods are evaluated in terms of a sensitivity factor which constrains 
the maximum white noise array gain. Results show that both methods result in sidelobe 
levels that range from 20 to 35 dB lower compared to a conventional beam pattern with 
uniform aperture weighting and that the design procedure is applicable to beam patterns 
steered to both true broadside and to off-broadside directions. 


Thesis Supervisor: Dr. Arthur B. Baggeroer 
Professor of Ocean Engineering and Electrical Engineering 
Massachusetts Institute of Technology 





Acknowledgements 


My two years as a graduate student have been an enjoyable and rewarding period. 
Many people have aided me in the course of my studies, more than 1 can reasonably 
mention here. However, there are a few people who deserve special mention. 

First and foremost, 1 must thank my advisor, Dr. Arthur B. Baggeroer. Without 
his invaluable guidance, encouragement and patience, I would not have been able to 
successfully complete my studies. I must also thank Art for his concern and interest in 
my career. 

A special thanks also goes to Eddie Scheer, computer wizard and musical guru. Eddie 
had to contend with numerous computer related questions and requests while maintaining, 
for the most part, his jovial attitude. In particular, Eddie was an available outlet when a 
break was needed and was an immense help when he typed a portion of this paper. 

I am also indebted to Rob Fricke. Rob was the one person who was always available 
and always willing to discuss my studies and unrelated topics. 

I must also recognize my fellow Joint Program adventurers. Although many of us had 
(and still have) the same goals, we have achieved those goals through different avenues 
and in doing so, have acquired a shared base of knowledge not found anywhere else. 

Iam most grateful to the United States Navy in conjunction with the Oceanographer 
of the Navy for allowing me this valuable opportunity and for provided the financial 
backing. 

And last, I must recognize the love and support I have received from my family. My 
wife, Laurie, and my daughters, Brittany and Brianna, have kept me in touch with reality 
and my true priorities in life. 





Contents 


1 Introduction 9 
MS aero tn dE a e e a Denn لځعهکښ‎ A 9 
EN N ee AOL 10 
NN ی‎ eee era 11 

2 Problem Development 13 
|1 ٢1 171 ۱2011 کال له :1 هه هد‎ vas لمکا هه ره هه‎ es See a 13 
DENE vray Processing Approach .........2.2 m mn nen 15 

n ena Node o a E r a a P a a DOD ca Un ol هه‎ 15 
Dx spectral Covariance Matrix . 00... oe 17 
E S A A el ee 21 
A A e اوو سوه ور‎ BOD ORS b 24 
DAS Spatial Scanning aa ee a T 25 
ew Beam Pattern Design Method . . ............<.......<... 29 
A RAS ee ee ee 34 
Pee N Te TO e a m ai me 35 
2.4 BeamformerPerformanceMeasures....................... 35 
e enva lues a a a A İİ e e e 35 
د د‎ SESI ٢٢٥ C.L CCE SN LTD 37 
O E و‎ e rcs me EE oen ee oe Oh a 38 

3 Beamforming Results 40 

re en ee Y 40 





Pe en es 
3.2.1 Method A 
3.2.2 Method B 


3.3 Steered Results 


pa Se مو‎ e e ۵ O e e ۷ 3^ 3.» «^ 6۰ 6 à DDE OE O ONE € . اه‎ 
هده د يد‎ o òo ē òo òo ?*? 9? e... ي‎ a e. o.  @  « 


3.3.1 Method A 
3.3.2 Method B 


e » ò> ò> ò òo 8 ē òo ē GE ^» ē ò% ē ò> ē ò> ^? ē ò> ē ò> ē ò> Gg ^? ها‎ e. هوت ها‎ 


€ > ò% ò e e e» e...» e ò ò ò> ò e» ò> ده ھت ھت هده هت‎ ò> ò ò% ò% ò% ò% e» € 


4 Summary and Conclusions 
A Spectral Covariance Matrix Inversion 


B Dolph-Chebyshev Weights 


T4 


T8 


81 





List of Figures 


2-1 
2-2 
2-3 
2-4 
2-5 
2-6 


2-7 


2-8 
2-9 


2-10 Two Dimensional Cross Section of Penalty Function 


3-1 
3-2 
3-3 
3-4 
3-0 
3-6 


TTA EO ECEN e Ne ee e 
Array Coordinate System and Dimensions + os ی مه‎ a a 
Element Numbering Scheme o. ue vw ie A Ee s 
Specural@ovanlance Wiathixs «cits Locos. ve e ee o I E E LE RUM uus 
System Architecture for a 200 Element Array . .. 2-2 e 58 5150 
Projected Scanning Vector Locations. Every other scanning increment in 
elevation and azimuth is shown. The shading illustrates the variability in 
azimuth of the spatial separation between scan directions relative to the 
spatial separation at broadside resulting from the uniform angular scanning 
EA ita ers ههاو هره مق غار نو‎ oa 
Null Placement Effects. For a source (solid line) incident at some k, the 
MVDP beamformer optimally places a null in the array response pattern 
(dashed line) in the direction of the source incident at k. 


Three Dimensional Mesh Plot of Penalty Function 


. . . e e Gg 4G 8 8 e O o . 


Overhead View of Penalty Function 


€ @ o% Gg 9 ò% G 4G % G .—.—.. gg a ^*^ ^» G o ə% 


و o. e e‏ ما ها o‏ ما ها ها @ . 


array output for plane wave case 
conventional beam pattern steered to O° elevation and 90° azimuth .. . . 
conventional beam pattern steered to 0° elevation and 102° azimuth 

Eigenvalue Performance Results for Case 11A 


. . oe ها ها ها ها ها‎ o 8... . ها‎ e 


array output for case 11A, Amin = o~ 


SA è% s o ọọ y o a e e 9 S « e e G o 4 


beam pattern for case 11A, Amin = 107° 


u. 8 90 e 4G 6 . e. ii 4. @ e > s.s. ð: s € 


14 
14 
19 
20 
23 


43 
44 
47 
48 
49 





3-7 square magnitude of element weights, case 11A, Amin = BO SERIE. 
3-8 corrected phase of element weights, case 11A, Amin = 11 a NN 
3-9 Eigenvalue Performance Results for Case 9B................. 
3-10 array output, case 9B, Amin = A LASER O. 
3-11 beam pattern, case 9B, Amin = 1 NS a ey a eret 
3-12 square magnitude of element weights, case 9B, Amin = LORE rM 
3-13 corrected phase of element weights, case 9B, Amin = NOs eb, are, 
Su array output, case 5B, Amin = 1O0 * 00. um 
3-15 beam pattern, case 5B, Amin = E e EA at te cs ANGE AĞ 
3-16 square magnitude of element weights, case 5B, Amin = 107% ........ 
3-17 corrected phase of element weights, case 5B, Aj, — 107 *$ ......... 
3-18 array output, steered method A penalty function .............. 
3-19 beam pattern, steered method A penalty function. ............. 
3-20 square magnitude of element weights, steered method A penalty function 

3-21 corrected phase of element weights, steered method A penalty function . . 
3-22 array output, steered method B penalty function .............. 
3-23 beam pattern, steered method B penalty function. ............. 
3-24 square magnitude of element weights, steered method B penalty function 


3-25 corrected phase of element weights, steered method B penalty function . . 


4-1 Plot of First Sidelobe versus Center Region Width, Method A. Conven- 
tional main lobe widths in terms of angular separation between first nulls 
are approximately 24° in azimuth and 50° in elevation... ........2.. 

4-2 Plot of First Sidelobe versus Center Region Width, Method B. Conventional 
main lobe widths in terms of angular separation between first nulls are 


approximately 24° in azimuth and 50° in elevation 


e 9 ò% e ò e e o e e ee ò ò o 


هه  .‏ ما ها ها ها »^ * s,‏ . . 


B-1 Dolph-Chebyshev square magnitude element weights 


61 
62 


65 


67 
68 
70 
71 
12 
73 





List of Tables 


2.1 


3.1 


3.2 


3.3 
3.4 


3.0 


A.l 


Summary of Computational Requirements for Different Scan Patterns .. 27 


Penalty Punction Parameters 2.2.0 375 24x: a NN لاله‎ 45 
Minimum Eigenvalues and Associated Dynamic Ranges Tested for Method 

A tcn E e a 46 
Method A Case Study Summary 244309933339 9 RO Rmo 47 
Minimum Eigenvalues and Associated Dynamic Ranges Tested for Method 

|| ار‎ 11 1101015 a San as Soe See OUR US we es 53 
Method B Case Study Summary . >... se se sos so s soo ss ese. 54 
Computational Times for Matrix Solution Methods ............. 79 





Chapter 1 


Introduction 


1.1 Background 


The problem of beamforming with multidimensional arrays is a very important prob- 
lem with numerous applications, yet it is not well understood. Whereas one dimensional, 
equally spaced arrays have many formulations and solutions based on extensive temporal 
analysis, arrays with added dimensions and nonuniform spacing introduce unique prob- 
lems and possibilities which have no direct temporal eguivalent.|1| In terms of the bearing 
estimation problem which is addressed in this investigation conventional beamforming 
can be described as the shading or windowing of array sensor data to obtain a desirable 
directional scan pattern, or beam pattern. One dimensional windows are numerous and 
well documented.[2] In two dimensional and higher dimensional cases, fewer windows are 
available.[3] The majority of the multidimensional windows are based either on a perfect 
circular based symmetric geometry (a perfect sphere or cylinder) or on the extension of 
one dimensional windows to higher dimensions which routinely requires uniformly sam- 
pled apertures.[4,5] Unfortunately, all multidimensional arrays are not completely circular 
and/or do not have regular sensor spacing. The impact of array geometry, particularly in 
three dimensions, on the spectral estimation and bearing estimation problems is not well 
understood.[6] Therefore, the purpose of this investigation is to develop a beam pattern 
design procedure that can be used on general multidimensional irregular arrays which 


includes the effects of the array geometry. The design procedure uses a data-adaptive, 





minimum variance method that can be applied to arrays with nonuniform spacing. The 
minimum variance procedure is Capon's High Resolution Frequency Wavenumber Esti- 
mation Method more commonly known as the Maximum Likelihood Method (MLM).|7| 
The test array, which has proposed applications in an undersea environment, consists of 
a three dimensional sonar array located on an arc of a cylindrical surface. The primary 
objective is to design beam patterns by processing array data with the MLM beamformer. 
The source of array data is the design tool and is contained in a special spectral covariance 
matrix termed a “penalty function.” Generation of the penalty function is the important 


issue of the design procedure and is the focus of this investigation. 


1.2 Objectives 


As previously stated, the primary objective is to develop a beam pattern design pro- 
cedure using an MLM beamformer that can be used on irregular arrays where traditional 


aperture shadings are not applicable. This primary objective results in two issues: 


e can the MLM beamformer be effectively used in a design role? 
e can the MLM beamformer control sidelobes? 


Both of these issues are examined using the penalty function concept. Selection of the 
proper penalty function results in a beam pattern with low sidelobes when compared 
to a conventional beam pattern generated from a uniform aperture shading. Different 
penalty functions result in beam patterns with different sidelobes. The selection of dif- 
ferent penalty functions gives the designer control over the procedure. In applying this 
procedure the designer obtains a qualitative perception of the array’s performance in 
terms of the array geometry and element spacing which aids in the design process. 

To aid in the analysis and the evaluation of the primary objective, quantitative mea- 
sures of array performance in terms familiar to conventional beamforming are included. 
The two performance measures relied upon are directivity index and sensitivity ratio. 
Sensitivity ratio is a relative measure of white noise array gain or array sensitivity. The 
quantitative measure used to describe the penalty function in terms of array performance 


is the minimum eigenvalue of the penalty function. One additional aid is a comparison of 


10 





the shading obtained by the design procedure with the aperture shading resulting from a 
Dolph-Chebyshev design for a uniform rectangular array. The quantitative performance 
measures and Dolph-Chebyshev comparison study enable an analysis and evaluation based 
on proven conventional methods and parameters. 

The array application used in this investigation deals with a large array composed of 
200 elements in a three dimensional sound field. Arrays with a large number of elements 
(as is the case here) lead to computationally complex and time intensive array processing 
requirements; therefore, a significant issue concerns obtaining an implementation which 
reduces or limits the computational complexity and processing times. In order to limit 
the processing requirements, the aspects of directional (baffled) hydrophones or sensors, 
array structural shading and the estimation process of the spectral covariance matrix are 
not considered. The array structure is transparent to sound from any direction and the 
spectral covariance matrix is assumed to be available. The procedure is based on a plane 
wave development. Since the MLM requires a matrix inverse, even the use of these sim- 
plifications still leads to a time intensive process. In an effort to further reduce processing 


times, the investigation highlights a comparison of two matrix solution methods: 
e Gaussian elimination 


e a Toeplitz bordering approach based on Levinson’s method.|8| 


1.3 Overview 


The remainder of this paper is divided into three chapters and two appendices. 


Chapter 2 describes the problem and is broken into four additional sections. 
e Section 2.1 describes the array geometry. 


e Section 2.2 defines the signal model, develops the spectral covariance matrix and 
Its properties, reviews conventional beamforming and MLM beamforming, defines 


array gain and develops the spatial scanning problem. 


e Section 2.3 develops the penalty function. 
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e Section 2.4 develops and defines the beamformer performance measures. 


Chapter 3 contains the beamforming design results. The design results are split into 


three sections. 


e Section 3.1 includes beamformer verification results which are useful as points of 


comparison. 
e Section 3.2 contains extensive results on broadside beam pattern designs. 
e Section 3.3 contains results of steered beam pattern designs. 


Chapter 4 contains a summary. 

Appendix A contains timing data resulting from the investigation of the Gaussian 
elimination and Toeplitz bordering matrix solution methods. 

Finally, appendix B contains a comparison of the aperture shading design results with 


the traditional Dolph-Chebyshev design procedure. 
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Chapter 2 


Problem Development 


2.1 Array Geometry 


The array used in this investigation has proposed applications on submersible bodies. 
It is located on a 60 degree arc of a cylindrical surface. The array is shown in figure 2-1 
and is to be operationally oriented as depicted. The array has intrinsic structure but is 
neither linear nor circular. The array is to be applied in a manner which reduces the 
problem to processing only the half space on the same side of the cylinder on which the 
array is located. 

The array consists of 200 elements arranged in ten rows (or lines) of twenty elements 
each. Considered separately, each line of twenty elements is identical. The array is 
symmetric about the x and z axes for the coordinate system defined in figure 2-2. This 
structure leads to some important results when considering the spectral covariance matrix. 

The design wavelength is defined as the wavelength (A) that determines the element 
spacing. All of the work performed in the investigation of this array is accomplished at the 
design wavelength; therefore, design wavelength and A are interchangeable. The element 
spacing is specified to be one quarter of the design wavelength or A/4. All dimensions of 
the array (including the cylinder radius) are related in terms of the design wavelength. 

The structure of the array can be compared to a simpler two dimensional planar array. 
This comparison is valuable because uniform planar arrays are well documented. In order 


to prove this comparison is justified to first order, the extent the actual three dimensional 
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Figure 2-1: Array Geometry 


6 = 0.2282 — L = 4.754 


Figure 2-2: Array Coordinate System and Dimensions 
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array deviates from a true uniform planar array must be determined. If the array is 
considered to be a structure of ten 20 element line arrays , then the length (L) of the line 
arrays is 4.754. If the angular separation of the line arrays 1s uniform and the arc length 
of this angular separation is approximately linear and equal to 1/4, then the relationship 
between R and A can be determined by equating the arc length of the 60 degree sector to 


the arc length enscribed by the 10 line arrays. This results in 


7 A 
—R = 9— 21 
3 4 ( ) 
or 
pu p (2.2) 
An 


By applying some simple trigonometry, the planar height (H) of the array as projected 
onto the x-z plane ıs equal to the radıus R, or 2.15X. The maximum deflection of the array 
(or a measure of the curvature from the x-z plane) is approximately equal to 0.2881. The 
projected vertical spacing 1n the array center is equal to 0.25A and is approximately equal 
to 0.2231 at the array top or bottom. el explanation of these dimensions is 
ıllustrated in figure 2-2. The projected array has uniform spacing along the axial (or z) 
axis and nonuniform spacing in the x direction (although the spacing is nearly uniform.) 
Therefore, the three dimensional array can be approximated by a uniform rectangular 
array if the effects of curvature are neglected. The rectangular array approximation 
is better in the axial direction; therefore, the expectation is that design results in the 
axial direction will more closely resemble results of the uniform rectangular array. The 
rectangular array comparison provides a reasonable model which is useful in establishing 


an intuitive feel for the expected performance of the test array. 


2.2 Array Processing Approach 


2.2.1 Signal Model 


One purpose of an array and its associated processor is to provide an estimate of 


the directional spectrum of a space/time stochastic process. The directional spectrum 
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information is contained in the freguency wavenumber function. To understand the ar- 
ray processing system the signal model must first be defined. The space/time stochastic 
processes being considered are assumed to be temporally wide sense stationary (defined 
in reference [9]) and spatially homogeneous (defined in reference [1].) The wide sense 
stationarity assumption implies that disjoint frequency bands are uncorrelated. The spa- 
tial homogeneity assumption implies that disjoint wavenumber bands or directions are 
uncorrelated, i. e. no correlated multipaths. Both the wide sense stationarity and spatial 
homogeneity assumptions are not universally valid in the ocean environment, but can be 
closely approximated in practice by proper choices of time and space scales. Although 
time domain analysis is not addressed in this study and the actual physical dimensions of 
the array are not known, the wide sense stationarity and spatial homogeneity assumptions 
are presumed to be valid based on anticipated time and spatial scales. The space/time 


stochastic processes are represented in a Stieltjes integral form [10] as 


where t is the time argument, z is the vector representing the spatial position and k is 


the three dimensional wavenumber vector.! 


For the wide sense stationary and spatially 
homogeneous processes considered, X(f,k) defines a Fourier transform representation for 
a plane wave of frequency f and wavenumber vector k. Thus z(t,z) is composed of a 


superposition of plane waves X(f,k).{1| For the three dimensional ocean environment, f 


and k must satisfy the dispersion relation 
۱ nie (2.4) 


Equation 2.4 results from the three dimensional wave equation where c is the speed of 
sound in water. This dispersion relation corresponds to free space propagation. The 
space/time covariance function K(t,t—7,z,z— Az) is the cross correlation of the stochastic 
process at times t and t — 7 and spatial positions z and z— Az. All processes are assumed 

‘Lower case letters are used to depict temporal representations while upper case letters depict transform 


representations. Similarly, upper case bold face letters represent matrices and lower case bold face letters 
represent vectors. 
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to be zero mean stochastic processes. HU t —7,z,z— Az) is determined by 
K,(t,t—1,z,2— Az) = Elzx(t,2)3*(t-7,2- Az)) 25) 


where * denotes complex conjugate and E represents the ensemble expectation operator. 
The frequency domain representation of the space/time covariance function is the spectral 


covariance function S,(f,z,z— Az) which is defined as 
S,(f,z,z — Az) = J K.lt,t-r,2,2- Az)e ’?"!"dr (2.6) 


S,(f,z,z — Az) represents the cross spectra between positions (corresponding to array 
elements) z and z — Az. The frequency wavenumber representation of the space/time 


stochastic process, called the frequency wavenumber function P( f, k), is defined as 


P(f,k) = J J [ > S.(f,2,2— Az)ekAZdAz) (2.7) 


Since f and |k| are related by equation 2.4, fixing f reduces the problem to determining 
the direction of propagation of |k|. This is the approach used throughout this study; f 
is fixed and then the half space in which the array resides is scanned to determine the 
direction of propagation parameterized in terms of an bearing in azimuth and an elevation 
angle. The estimate of the direction of propagation or directional spectrum, embedded in 


the function P(f,k), is the desired output of the array processor. 


2.2.2 Spectral Covariance Matrix 


An examination of the spectral covariance matrix is included at this point because it 
alds in the understanding of the beamforming procedure and also reveals some pertinent 
characteristics of the matrix operations which are encountered. 

The spectral covariance matrix can be decomposed into two components, one coherent 


and one incoherent, by setting (dropping arguments for brevity)? 
m 1 2 
S — osdd ap INQ (2.8) 


24 represents complex conjugate transpose 
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where 
el Ke "Zi 


d= ۱ )2.9( 
elke Zn 

is the geometric phase or direction vector of a plane wave signal propagating with wavenum- 
ber vector k, and signal strength 03. Q represents the normalized (trace Q =N) noise 
cross spectral density matrix of element outputs with noise strength o. If the noise 

consists solely of uncorrelated components, then Q reduces to the identity matrix I. 
The coherent input places the restrictions on f and kg found in equation 2.4. The 
coherent input used in this study consists of plane waves which are modeled at time t and 


position Z as 


x(t, 2) — a(f)e(2r11-&e2) (2.10) 


where a(f) is a frequency dependent amplitude function. The coherent component of the 
spectral covariance matrix resulting from the plane wave input (determined by applying 
equations 2.5 and 2.6) is 

oz ddl = ake Tee Otis Vij (2.11) 


where g% is frequency dependent and Az; ; = z; — z; corresponds to the vector difference 
between the f and J elements. The exponential argument k, - Az; ; represents the 
phase delay due to propagation of the plane wave between the #** and j* elements. 

The first observation deals with the structure of a matrix formed by equation 2.8. 
Matrices of this form are Hermitian symmetric or simply Hermitian. A Hermitian matrix 
is defined by S = St. The second observation results after examining equation 2.11. 
The coherent input is dependent only on the difference in element positions. This spatial 
homogeneity leads to an embedded Toeplitz structure. A Toeplitz matrix is defined as a 
matrix whose ijth element s(t, 7) is a function of (i-j) and thus has identical elements along 
the main diagonal and each of the subdiagonals. (Reference [11] contains an excellent 
summary of Hermitian and Toeplitz matrix properties.) To fully see the structure of 
the matrix formed by equation 2.11 for the test array, consider an element identification 
scheme where the elements are numbered consecutively along each successive row (or line) 


starting at the array top. This numbering scheme is shown in figure 2-3. The spectral 
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Figure 2-3: Element Numbering Scheme 


covariance matrix resulting from this numbering scheme is shown in figure 2-4 where the 
matrix is partitioned to indicate specific sub-matrices. Each s, corresponds to the cross 
spectral covariance of the kth element with the [th element. Each sub-matrix corresponds 
to the spectral covariance formed between elements located in two of the line arrays. The 
main diagonal of sub-matrices is the result of forming the spectral covariance of each line 
array with itself; the main diagonal sub-matrices are all identical because each line array 
is identical and the coherent input is perfectly correlated. 
S. can also be represented as a two level matrix where to the first order level 
S11 Siz ° S140 
ao ۱ )2.12( 
$101 $1010 
Each sub-matrix S;; is the cross covariance between the ith and jth line arrays and 
is defined in figure 2-4 by each partition. For example S;2 is the sub-matrix in the 
partitioned first row and second column of figure 2-4 and is formed from the line array 


with elements 1 through 20 and the line array with elements 21 through 40. The spectral 


covariance matrix has several distinctive characteristics. 


e S, is NxN (or 200x200 for the test array.) 
e Each sub-matrix, S;;, is Toeplitz and Hermitian 


e S,, to the first order level, is Hermitian but not Toeplitz because of the vertical 
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structure of the array. However, the main diagonal isa constant (to the first order 


level.) 


e The addition of the incoherent term will not change the embedded Toeplitz structure 


of the sub-matrices. 


Matrices with the form of S, are frequently classified as a two level matrices. Ref- 
erences {12] and [13] include a complete discussion on the classification of matrices in 
this manner. S, in this case is classified as a Hermitian-Toeplitz two level matrix, but 
by permuting corresponding rows and columns, the matrix can be transformed into a 
Toeplitz-Hermitian two level matrix where the first level is Toeplitz and the second level 
sub-matrices are Hermitian. This is important because there currently exist efficient and 
documented algorithms which can invert a Toeplitz-Hermitian two level matrix.[13] As is 
shown later the spectral covariance matrix needs to be inverted. Details of the matrix 
inversion selection process are included in Appendix A. (This selection process played a 
very important role in this investigation; the time intensive inversion process had a direct 


relationship to the processing time.) 


2.2.3 MLM Development 


In the context of this investigation array processing and beamforming are synony- 
mous. There is much literature on beamforming, both conventional and data adaptive. 
Complete reviews of processing methods are readily available.[14] Since the purpose of this 
study is to develop a beam pattern design procedure that can be applied to nonuniform 
multidimensional arrays, the choice of usable methods is limited to the MLM beamformer. 
The procedure can be more accurately described as a minimum variance, distortionless 
processor (sometimes referred to as MVDP).[6] From this point on, the beamformer will 
be referred to as the MVDP because it more accurately describes the processor. What 
follows is a development of conventional beamforming and the MVDP. 

w(f|kr) is the aperture weighting function with the processor steered to kr.” For the 


sampled aperture consisting of an N element array, the it? component of w is the complex 


“The development is the same as used in reference [1]. 
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conjugate of the 1** element's weighting (w;). 


w(f | kr) — | wi(f|kr) (2.13) 


The aperture weighting function, or weighting vector, can also be considered to be an n- 
dimensional taper or shading function for an n-dimensional array. It should be observed 
that w is dependent on the steer direction and can be complex. The scanning vector 


e(kr) steered to kr is defined as 


giEr ti 


e(kr) — 3 (2.14) 


elkr-tn 


where again z; is the 1% element’s position vector of an N element array. The scanning 
vector can be interpreted as the plane wave phase delays which must be applied to the 
array elements to steer the array in the direction of kr. The scanning vector also has an 
interpretation of being a replica field or the field at the array produced by a source at 
some position in the water column. 

The output of the array beam steered to kry is defined in conventional beamforming 


as 


N 
Y (f,kr) =wX => w;(f|kr)X(f, z:) (2.15) 


(=1 
This system architecture is depicted in figure 2-5. The array response pattern W(f, k[k7) 


is given by 


N 
W(f,k|kr) 2 w'e 2 Y w(f|kr)e* *: (2.16) 


=1 
and represents the response of the array (steered in the direction of kr) over the aperture 
to a plane wave signal. If k is fixed and the function (f, k|kz) is evaluated over a spatial 


region, then the beam pattern B is obtained where 
B(f,0,9|kr) — W(f, k|kr)|— 2۳ a(0,0) (2.17) 
where A is the wavelength and a(0,$) is a unit vector specified in a spherical geometry. 
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Figure 2-5: System Architecture for a 200 Element Array 


Because P(f,k) is a measure of power, this leads to a representation of the response of 
the array as a function of two array sensor locations. This will result in a second moment 
representation for the array output and requires a second moment input. The array input 
is the spectral covariance function S;(f,zi,z;) as defined by equation 2.6. The output of 


the beam steered to kr is the power density spectrum S,(f|kr) and is given by 


S,(f|kr) 7 w'(f|kr)Sz(f, zi z;)w(f|kr) (2.18) 


Sy(f|kr), as defined by equation 2.18, is the output of a conventional beamformer and 
is the estimate of P(f,k) with the array steered to kr. The MVDP is obtained by 
minimizing the response of equation 2.18 (minimum variance) subject to unity gain in the 


look direction (distortionless response.) The MVDP solutions are [1,7,15] 


wMvpP(f|kr) — Fie (2.19) 
el -le 

Wmvpr(f, k|kr) = A (2.20) 

Sy,mv pP(f|kr) = E (2.21) 


23 





The MVDP method adaptively finds different aperture weightings for different input signal 
fields such that sidelobes are adjusted to reject interfering coherent noise sources and 
the incoherent noise. This method is optimum in the sense that it minimizes the array 
response over the aperture of the noise field (the input signal not in the beam) by design 
and application of a spatial filter. This spatial filter estimates the power in a plane wave 
propagating with frequency f and wavenumber vector ky in the presence of interfering 
plane waves and sensor noise.{16] To fully define the directional spectrum, the signal field 
must be scanned over the appropriate angular space (which, for this study, is the half space 
in which the test array is located.) If only uncorrelated noise is present, then uniform 


weighting with the conventional beamformer is optimal. 


2.2.4 Array Gain 


Array gain is the improvement in signal to noise ratio due to beamforming. This is 
a result of finite beam widths rejecting noise and signals not in the beam. Ideally, for the 
signal of interest, signal rejection is small while the noise rejection is large resulting in a 
large array gain. The array gain as such is defined as 


[wi dl? 
= 2.22 

w!Qw (222) 
which is simply the ratio of the signal response to the noise response. In the context of 


the MVDP processor, the optimum array gain is [1,17] 
Gor < CO le (2.23) 


For the case of spatially white or sensor noise with omnidirectional sensors, the array gain 


becomes the white noise gain or 


$412 
MS ای‎ <N (2.24) 


where N is the number of sensors. Equality occurs when the processor is perfectly matched 
and there is uniform weighting. The processor is defined to be perfectly matched to the 
signal directional characteristics when e = d where d is defined by equation 2.9.[17] 


Mismatch occurs whenever e 4 d and may, for example, be caused by non-plane wave 
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propagation, errors associated with the sensors, sampling, quantization, and steering off 
of the signal direction. When the errors due to mismatch are uncorrelated, the sensitivity 


(S,,) of the array gain due to signal mismatch is related to G,, by [18] 


1 
e > 2.25 
ae 2) 


Used in this context, Sy or Gy is a measure of robustness. Using wy pp calculated 
in equation 2.19 for a perfectly matched processor with S, = I, G, equals N and the 
beamformer is robust. (wav pp for this case is a normalized, uniform weighting.) Any 
deviation from N indicates a degradation in beamformer robustness. 

One additional expression for the array gain which turns out to be useful relates G in 


terms of eigenvalues and eigenvectors of Q. The expression is 





(2.26) 


N 2 
İdiv;! 

G= > a 
i=] : 


where A; and v; are the eigenvalues and eigenvectors, respectively.|18] This expression 
shows that G is dominated by small ratios of |d!v;|?/A;. This effect is not well understood 


because when X; is small, dtv; is usually small also. 


2.2.5 Spatial Scanning 


The remaining feature to be specified is the spatial scanning problem. The spa- 
tial field, as previously defined, is the half space in which the array is located. To fully 
define the directional spectrum, the spatial field must be scanned or sampled. Equa- 
tion 2.14 defines the method yet does not answer the question, “how finely must this field . 
be sampled?” The spatial sampling interval is defined in terms of Nyquist’s sampling 
requirements. This requires sampling at one half the natural beam width of the array 
which is equivalent to sampling at twice the highest frequency of the field or spatially at 


A/2. Three factors must be considered when determining the sampling requirements. 


1. As described in reference [19], finer sampling results in a visual display with signifi- 
cantly more content. This additional content is valuable when evaluating hard copy 


or video outputs of the directional spectrum. 
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2. If the natural beam width is determined by assuming the array 1s approximately 


rectangular, then error will be introduced into the values obtained. 


3. The MVDP beamformer is susceptible to serious signal suppression effects due to 
signals arriving form directions between the processor beams. Minimizing this effect 


requires more closely spaced scan directions.|17] 


Consideration of these three factors leads to the selection of a scanning increment of 
1/4 the natural beam width (BW). Verifying this is the proper decision requires the 
evaluation of 1/2 and 1/8 BW scanning increments. There is one additional consideration 
when determining the proper scanning increment; it involves the computational burden. 
To illustrate this, the natural beam width of the array must be determined. 

For the test array, BW is resolved into two directions called elevation (also known 
as depression/elevation or D/E) and azimuth. These directions arise naturally when 
anticipating an operational scenario; elevation is an angular measure off of a horizontal 
plane (+90 degrees to -90 degrees) and azimuth is an angular measure fore and aft of a 
vertical plane through the center of the array (0 degrees forward to 180 degrees aft). To 
determine BW the assumption that the array is approximately uniform and rectangular 
is applied. BW for a uniform rectangular array in the x direction is [20] 


BW, = — 
E (2.27) 


For the azimuthal direction, Lazimuth = 4.75A and 
BW azimuth = 0.0677 x 12° (2.28) 


This leads to an azimuthal scanning increment of 3 degrees. In elevation (where the 
rectangular array approximation is least valid), Letevation = 2.154 and 
o 
BWV elevation = SN 265 (2.29) 
This leads to a scanning increment in elevation of 62 degrees. 
The scanning increment is implemented using uniform angular separation. The scan- 


ning system is best described using a global coordinate system of latitude and longitude 
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54 | 118 | 6372 
1593 
420 






Table 2.1: Summary of Computational Requirements for Different Scan Patterns 


where latitude corresponds to azimuth and longitude corresponds to elevation. The pro- 
cedure implementing the spatial sampling is to choose La latitudes separated by the az- 
imuthal scanning increment. Then along each latitude, L. longitudes are chosen separated 


by the elevation scanning increment. The total number of scans, L, is 
CSE: ba (2.30) 


Steers in azimuth (latitude) to 0° and 180° are independent of the elevation angle. To 
eliminate duplicate scans and place a scanning direction at 90° azimuth, a 3° offset is 
incorporated at 0° and 180? azimuth. This pattern results in Lg = 59. To place a 
scanning direction at 0°, a 31 offset is incorporated at +90° and —90° elevation; this 
results in L, — 27. Therefore, the total number of scans is L = 1593. The scanning 


pattern as projected onto a plane surface is shown in figure 2-6. 


Now consider the computational requirements of this scanning pattern. To define 
the directional spectrum equation 2.21 must be calculated L times, once for each scan 
direction. For an N element array, equation 2.21 represents an N by N system of equations. 
Solving an N by N system of equations L times is a time intensive process. Likewise, to . 
determine the beam pattern, equation 2.20 must be calculated L times. Fortunately, the 
N by N system of equations only needs to be solved once for the array steer direction. 
Table 2.1 lists the scanning requirements for 1/8 BW, 1/4 BW and 1/2 BW scanning 
patterns. Based on these computational requirements the desire is to sample the spatial 
field at the fewest number of points consistent with adequate beamformer performance. 

The scanning pattern used is not ideal. The algorithm implementing this pattern is 
very simple; with simple modifications 1/2 BW and 1/8 BW scanning patterns can be 


tested. However, the use of uniform angular scanning as implemented does not translate 
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Figure 2-6: Projected Scanning Vector Locations. Every other scanning increment in 
elevation and azimuth is shown. The shading illustrates the variability in azimuth of the 
spatial separation between scan directions relative to the spatial separation at broadside 
resulting from the uniform angular scanning implementation. 
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into uniform spatial scanning. Whereas the spatial scanning is uniform in azimuth, the 
scanning in elevation is not spatially uniform, particularly in regions close to 0% and 
180° in azimuth where the sampling in elevation is very dense. Figure 2-6 illustrates 
the differences in spatial scanning increment in elevation relative to the broadside spatial 
scanning increment. If this method of scanning is to be practically implemented, the 
optimum way is to rotate the pattern 90? and limit the elevation scanning region to the 
prime elevation angles of 4-60? where no structural shading occurs. This results in fairly 
uniform spatial scanning, but over a limited spatial region. 

In summarizing the spatial scanning problem, a simple method that suits the purposes 
of the idealized synthetic investigation of the test array 1s presented. The scanning pattern 
occurs in uniform angular increments of 3? with 59 steers in azimuth and 62 with 27 
steers in elevation for a total of 1593 steers. The spatial increment corresponds to 1/4 
BW scanning. The simple scanning method is easily modified to allow the use of different 


scanning increments, but does not provide uniform spatial sampling. 


2.3 Beam Pattern Design Method 


In attempting to design beam patterns, the objective is to reduce the sidelobe levels 
compared to sidelobe levels characteristic of traditional shading methods. The beam 
pattern design procedure is implemented by the use of a “penalty function,” so named 
because the function penalizes beam pattern lobes in regions of high directional intensity 
while enhancing the lobes in regions of low directional intensity. To understand this 
function, the actual processes undertaken by the MVDP beamformer must be understood 
in a qualitative manner. 

The MVDP beamformer designs optimum beam patterns by minimizing the energy in 
the beam subject to the constraint of unity gain in the scan direction. The beam energy 


is represented as 
S(flkr) = | P(f,k)|W|?dk (2.31) 


which is interpreted as the integral over all wavenumber space of the frequency wavenum- 


ber function weighted by the square magnitude of the array response pattern (or beam 
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directional intensity 


k 


Figure 2-7: Null Placement Effects. For a source (solid line) incident at some k, the 
MVDP beamformer optimally places a null in the array response pattern (dashed line) in 
the direction of the source incident at k. 


pattern). The MVDP beamformer designs array response patterns based on the frequency 
wavenumber function of the incident field in order to minimize the beam energy. If a dis- 
crete source not propagating in the look direction is incident (resulting in large P(f,k) 
for that direction), then |W|? is reduced to minimize the beam energy. The MVDP beam- 
former optimally places nulls in the directions of off-look direction sources.[15] This null 
placement effect is illustrated in figure 2-7. If diffused directional sources are located in a 
finite sector, then a lower array response occurs in the direction of the sector. However, 
the placement of deep nulls results in higher overall sidelobe levels which increases the 
effects of sensor noise, the other source term which contributes to total beam energy. 
Conversely, minimizing the effects of sensor noise by lowering overall sidelobe levels re- 
duces the effects of nulls and increases the beam energy from off-look direction sources. 
The effects due to discrete directional sources and sensor noise are compromised until the 
minimum beam energy is obtained. 

Since the desired effect is a reduction in sidelobe levels, the primary component of 
the penalty function consists of a field of diffuse directional sources of varying intensity. 
The field of directional sources is modeled similarly to isotropic noise which is directional 
noise incident uniformly from all directions with uniform intensity. The penalty function 
directional source component consists of a modified isotropic noise field, incident uniformly 
from all directions but with nonuniform intensity. In terms of equation 2.31, the penalty 


function (defined by P(f,k)) in the directions of beam pattern sidelobes has high intensity 
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while in the direction of the beam pattern main lobe, has low intensity relative to the high 
intensity region. A three dimensional mesh representation of the penalty function is shown 
in figure 2-8. In regions of high intensity, |W|? is reduced. In regions of low intensity, |W|" 
can be large because the weighted integrated response 1s small. If the depressed sector of 
low intensity directional sources is centered about the steer direction and encompasses an 
area equivalent to the main beam region, the main beam of the resulting array response 
pattern will be large (supported by the unity gain constraint in the look direction) relative 
to the sidelobe region which 1s suppressed due to the uniformly distributed high intensity 
sources located in that region. The net effect is the desired result: a beam pattern 
with suppressed sidelobes. Figure 2-9 shows an overhead view of the penalty function 
centered about the steer direction. (The penalty function is always centered about the 


steer direction.) As shown, the penalty function consists of three regions. 


1. The center region, defined by width in elevation (A0¿) and width in azimuth (A¢,), 


represents the low intensity region 


2. The transition region, defined by width in elevation (A0¿) and width in azimuth 


(Adı), serves as the interface between low and high intensity regions 
3. The surrounding background region of high intensity 


The issues pertinent to penalty function design are the shapes and widths of the center and 
transition regions, the amount of stabilization, and effects of the steer direction. Several 
of these parameters are discussed now; others are addressed in a later section. 

The method used to input the penalty function into the array processor must take the 
form of a spectral covariance matrix. Using the same concepts as found in equations 2.8 
and 2.11, the penalty function PF is defined as 

PE oI SN (2.32) 
where D is the sensor noise strength serving as a stabilization factor and Sy is the 
directional component representing a modified isotropic noise field in terms of a sum 


of plane waves. The plane waves model this directional field at one half the scanning 
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Figure 2-8: Three Dimensional Mesh Plot of Penalty Function 
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Figure 2-9: Overhead View of Penalty Function 
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Figure 2-10: Two Dimensional Cross Section of Penalty Function 


33 





increment. Sy 1s defined as 


200 200 
=E ISU) Y Y Ener (2.33) 
¢=0 6=0 i=1j=1 


where S(f,0,¢) is a spatial weighting function acting on the plane waves and a,(0, $) is 
a unit directional vector defining the direction of propagation. I 4-0 29-0 Tepresents the 
sum over the scanning space consisting of the half space in which the array is located. Sy 
is simply the sum of plane waves spatially weighted to form the distinct sectors shown 
in figure 2-9. A two dimensional cross section of S(f,9,¢) is shown in figure 2-10. This 
illustrates the regions of high and low intensity. Two methods implementing the penalty 


function are now examined. 


2.3.1 Method A 


Method A uses the direct application of null placement to reduce beam pattern 
sidelobes. The high intensity directional noise field is formed with unit intensity while 
the low intensity region is formed with zero intensity (or no directional field). The sensor 
noise is then used to stabilize the process. The resulting directional noise field causes 
the sidelobes to be suppressed while the main lobe level is maintained by the unity gain 
constraint. 

Referring to figure 2-10, method A is implemented by setting Smaz(f,9,¢) = 1 and 
Ap = 1. In terms of equation 2.33, this requires 


E ie 0, $) =0 (2.34) 
S sito» 0, $) = erp|— e p|- E E د‎ (2.35) 
li, 0, $) =1 (2.36) 


This produces flat center and background regions interfaced by a Gaussian shaped tran- 
sition region. ło and ġo are the directions representing the appropriate center/transition 
region boundary. For a fixed number of increments in the transition region, the region 
dimensions and roll off’s can be controlled by adjusting the variance (o2 for elevation and 


g4 for azimuth) of the Gaussian function. 
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2.3.2 Method B 


Method B uses the opposite of null placement to reduce sidelobes. Instead of placing 
nulls in the sidelobe region, anti-nulls are placed in the main lobe region. This is accom- 
plished by forming the high intensity region with zero intensity (no directional sources) 
while forming the low intensity region with negative intensity sources. The low intensity 
region still has lower intensity relative to the high intensity region. The effect of sidelobe 
suppression is indirectly accomplished with method B. The negative intensity (or low in- 
tensity) region can support very large main beam levels since the integrated beam energy 
from equation 2.31 is small. This effectively enhances the main lobe. The sidelobes are 
suppressed when the unity gain constraint scales the large main lobe. The sensor noise is 
again used to stabilize the process; however, greater stabilization is required for method 
B than is required for method A. 

Referring to figure 2-10, method B is implemented by setting Smaz(f,6,¢) = 0 and 


A p — 1. In terms of equation 2.33, this requires 


Se 0, $) — —İ (2.37) 

ps 0, 2 Sn 2 
SB transition J,0,9) = - | لي مس‎ )2.38( 
Op takoota J 0, $) = 0 (2.39) 


The region shapes and dimensions are controlled in the same manner as in method A. 
In order to discuss the other issues, quantitative measures which enable an evaluation 
of those relevant issues must be determined. In particular, measures are needed to evaluate 


the effects of region widths and the amount of stabilization required. 


2.4 Beamformer Performance Measures 


2.4.1 Eigenvalues 


Equation 2.21 requires that some specific restrictions be placed on the spectral covari- 
ance matrix of the penalty function. The most significant requirement is that S, must be 


invertible. For a meaningful problem, this is ensured by requiring S4 be positive definite. 
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The approach taken to ensure this requirement is satisfied centers around the eigenvalues 
ol >... 

In section 2.2.2, it was determined that S, is Hermitian which results ın some special 
eigenvalue properties. All eigenvalues in this investigation are calculated using the Hermi- 
tian matrix routines available in the EISPACK eigensystem package.[21] Since eigenvalues 
of Hermitian systems are always real [22], a simple characterization of S; is obtained. This 
simple characterization has physical interpretation which is now pursued. 


Consider the eigenvalue problem 
Szvi — Aivi (2.40) 


where A; is the i? eigenvalue and v; is the eigenvector corresponding to A;. S, can be 
characterized as a measure of the energy of the received signals (similar to the eigenvalue 
interpretation of array gain.) This leads to an energy interpretation of eigenvalues where A; 
is a measure of the energy projected in the direction of v;. Small A; represent eigenvectors 
with low intensity. Low energy levels are associated with eigenvectors corresponding to the 
low intensity region of the penalty function. Again, this 1s consistent with the eigenvalue 
representation of array gain where small eigenvalues correspond to components with less 


noise.[17] Since S; must be positive definite, the smallest eigenvalue, Amin, must satisfy 
A >0 (2.41) 


Intuitively, a penalty function with finite (and positive) energy throughout the scanning 
space as measured at the sensor is specified since negative or zero energy systems do not 
possess physical interpretations in the ocean environment. 

With this interpretation for S,, positive definite matrices can easily be generated. The 
value of Amin 1s controlled by the stability factor. By generating a fixed coherent field by 
equation 2.33, Amin is varied by adjusting D in equation 2.32. 

While the requirement of equation 2.41 is essential, Amin must also be considered 
jointly with Amaz, the maximum eigenvalue. In particular, the ratio Of Amaz/ Amin is very 
important. This ratio characterizes the dynamic range of S,. A large dynamic range 
places more demands on the beamformer than a small dynamic range; therefore, dynamic 


range is an important issue when considering beamformer performance. Dynamic range 
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is applied in the area of array sensitivity and superdirectivity and is considered in the 


next section. 


2.4.2 Sensitivity Ratio 


The robustness constraint used in the investigation incorporates the principles associ- 
ated with the white noise array gain. Sẹ, related to G, by equation 2.25, is a measure of 
sensitivity to uncorrelated tolerance errors.|17] To ensure robust beamformer performance 
a constraint is applied to S,,.{18] The robustness constraint is applied to prevent superdi- 
rective performance. Superdirective performance is undesirable because a superdirective 
array is very sensitive to sensor noise and uncorrelated mismatch errors. Superdirectivity 
is characterized by extremely low values of Gy. 

The actual constraint applied couples G,, to a parameter called the array sensitivity 
ratio y which is defined as 

n = Nw'w (2.42) 

Equation 2.42 is derived by substituting for wary pp (as defined in equation 2.19) into 

equation 2.24. After some matrix manipulation (since S4 is Hermitian, then Sz! is also 


Hermitian), the result is 


[SERIE 
G = 2.43 
w,MV DP elsi. (2.43) 
Calculating wtw yields 
tg-2 
fee ><. 2 44 
(etSzle)? Cee 
For the perfectly matched beamformer, 
wiw = E (2.45) 
Gu, MV DP 


Comparing to equation 2.25, it follows that the magnitude squared of the element weight- 
ing vector is a measure of beamformer sensitivity (or robustness.) Since Gw is always less 
than or equal to N, equation 2.42 is the ratio of the uniformly weighted G,, to Gy Mv DP 
or simply a measure of how far the MVDP beamformer white noise gain deviates from 
N. If w represents a uniformly weighted aperture, then 7 equals 1 and the beamformer is 
robust. As the beamformer becomes less robust, y increases and the array becomes more 


sensitive. 
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The sensitivity ratio as defined in eguation 2.42 is a very convenient measure of beam- 
former sensitivity since w is available. But y must be related to superdirective array 
performance. Experience shows that superdirectivity occurs at sensitivity ratios in the 


range of 2.5 to 3.0; this leads to the constraint 
n< 25 (2.46) 


The value of y is controlled by varying the dynamic range of Sz. The addition of 0%, to the 
diagonal elements of S; adds 0%, to each eigenvalue and does not alter the eigenvectors.[18] 
Therefore, the inherent structure of the spectral covariance matrix is not changed. If 
Amin << 1, then adjusting Ajj, does not significantly affect Amaz thus controlling the 
dynamic range. Large dynamic ranges lead to superdirective array performance; therefore, 


ny is controlled by limiting the dynamic range of Sz. 


2.4.3 Directivity Index 


Until this point, the analysis of the array performance has been limited to an ideal 
ocean; one that has a perfectly coherent signal in incoherent noise. To examine the effects 
in a more realistic ocean environment, the requirement of incoherent noise is relaxed and 
the case of isotropic nose is considered. 

Directivity index is defined as the array gain of a perfectly coherent signal in isotropic 
noise. Array gain can also be defined as the ratio, in decibel units, of the signal to noise 
ratio of an array to the signal to noise ratio of a single element. The directivity index is 
obtained by assuming a perfectly coherent signal which has cross correlation coefficients 


between pairs of array elements of 


(e), 1 Vij (2.47) 


and isotropic noise with cross correlation coefficients 


sin(42|z; — z;|) 


p AS 
( n)i 7 25| م2‎ B z;| 


(2.48) 


where z; is, as before, the position of the 1%” element. Isotropic noise is defined as having 


the noise power per unit solid angle constant in all directions. Reference [1] describes it as 
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the superposition of plane waves propagating from all directions with uniform statistical 
level. Isotropic noise is commonly proposed as a first order model for ambient sea noise. 
Under these assumptions, the directivity index (DI) can be calculated by 


NON un 


DTO (q l >> o 2 
> 25; wiwj(Pn)i;j 


(2.49) 


where w; is the weight of the 1% element.[23] DI for a uniform rectangular array with 
dimensions corresponding to the test array 1s approximately 18 dB.[20] 

Use of DI is an attempt to obtain a performance measure based on a “more” realistic 
ocean than previously considered. Since “real” ocean noise is directional in both the 
horizontal and vertical and ocean multipaths cause reduced signal coherence, the use of 
DI has practical limitations. However, DI is commonly used and is useful as a more 


realistic measure of beamformer performance in the real ocean. 
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Chapter 3 


Beamforming Results 


3.1 Preliminaries 


Three preliminary issues must be resolved before the analysis of beam pattern designs 
can proceed. To aid in beam pattern analysis, a visual output is necessary. The relevant 
outputs used in this study are contour plots of the spatially scanned output, the beam 
pattern and the square magnitude and phase of the element weights. The contour plots 
are generated from a uniform grid of points. For the array output and beam pattern 
contours, the grid is 27 by 59. For the element weight contours, the grid is 10 by 20. 
Since the element weights are complex, it is convenient to plot the square magnitude of 
the weights and a corrected phase. The square magnitude of the weights are scaled to a 
maximum value of 1 (or 0 dB). The corrected phase (®;) of the 1** element has the effect 


of the propagation delay across the array removed. This is calculated by 


®; = ¢; - Sa, Zi (3.1) 
where gö; is the element’s actual phase, ap is the scan direction, and z; is the element’s 
position vector. All four plots have axes corresponding to the angular position in the 
spatially scanned half space. For the magnitude and phase plots, these axes relate the 
element positions in terms of theìr location in the scanning field. 

The next issue is the verification of the beamforming algorithm. The test used to 


verify proper operation is a plane wave incident at 0° elevation and 90° azimuth (true 
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broadside case) in the presence of sensor noise only. For this test the expected results 
are a sensitivity ratio of 1 (i. e. Gu,Mv pP — N = 200) and a uniform weighting vector. 
The MVDP beamformer results are as expected with DI — 39.2. The array output for 
this plane wave case 1s shown in figure 3-1. Observe that the array output correctly 
detects the plane wave. Figure 3-1 is also useful in that it provides some insight into the 
azimuthal and elevational resolution of the array processor. The beam pattern resulting 
from the plane wave input is the conventional beam pattern steered to 0° elevation and 90° 
azimuth; this is shown in figure 3-2. In figure 3-2, observe that there is unity gain in the 
scan direction as required in the MVDP derivation. Also, observe the first sidelobes (in 
both elevation and azimuth) are at -13 dB which is as expected for a uniformly weighted 
rectangular array. The second sidelobes in azimuth at -18 dB also correlate well with the 
rectangular array approximation. These two plots are useful as points of comparison for 
the broadside results. Results are also included for steered patterns, in particular, patterns 
steered to 0° elevation and 102° azimuth. For comparison purposes in anticipation of the 
steered results, the conventional beam pattern steered to 0° elevation and 102° azimuth 
is included in figure 3-3. 

The last issue concerns the sufficiency of 1/4 BW scanning. When a simple penalty 
function, set up to ensure directions between scanning increments are present, is input to 
the beamformer for the cases of 1/8 BW, 1/4 BW and 1/2 BW scanning, no observable 
effects due to incident plane waves between scan directions are present. Visually, the 
finer the scanning increment, the better (qualitatively) the output appears which is as 


expected. These results confirm the sufficiency of ¿BW scanning. 


3.2 Broadside Results 


The analysis of this investigation centers on a systematic evaluation of the effects of 
different penalty functions on beamformer performance (as measured by the sensitivity 
ratio and directivity index) and beam pattern structure. The parameters specifying the 
penalty function are listed in table 3.1 and are identified in chapter 2.3. The center region 


widths are easily defined; they are simply the widths of the plateau. The transition region 
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Figure 3-1: array output for plane wave case 


A 
Done 
. 











scale in dB 


WELL 
7 7 


Y 
i 
7 


trit 


95 A A me S 
VIS- - S 
و‎ - d 
8'er- - 

Zor c 

S'9e- - 

62€- - 

262 - 

952 - 

612 - 

Du 

gv - 

01۲ ~- 

£^ - 

teo ~- 

Le 





(seaiBap) uolyeAgje 


60 180 
imuth 


*. 


1 


0 
d 90° az 


14 
10n an 


e 


20 


1 


) 


80 100 
azimuth (degrees 
43 


60 


40 


20 
conventional beam pattern steered to 0° elevat 


Figure 3-2 





scale in dB 


Ss 
V] 


۷" FR 


70 ٢ 


€ €9- Mo|9q 
poss CES 
6vSg- - | 6S- 
Zug G US ; = ; 
v'9p- - LOS BR e مر‎ ATIII.: : وښ نو رب ۱۳ : په‎ 1 
CO ^ vQ9v- ی بت وس‎ 
oge- ° Zt E 7 aI ی‎ : : : 
866 - OSE- 

9'6c- BEG 
$56: " 96c- 
Lic © Sc 


691 PEE le: 
۸61۸ ” 69l- 
VG AS 
ev» و‎ 

ev 2۸005 





(see1Dop) uoneAeJe 


180 
imuth 


160 


140 
d 102? az 


lon an 


100 120 


th (degrees) 
44 


80 


azimu 


40 60 


20 
conventional beam pattern steered to O? elevat 


. 
* 


Figure 3-3 













center region width A0. 
ld > 
transition region width | elevation | A6; | 
A [a 


Table 3.1: Penalty Function Parameters 





widths are more difficult to define. The widths are controlled by varying o; and 05 1۳ 
equation 2.38. From this, A4; and Ağ; are defined in terms of the 3 dB down points 
along the major axes. All region widths are measured in degrees. The investigation 
concentrates on penalty functions centered about 0° elevation and 90° azimuth; the bulk 
of the results come from the study of the broadside case. This scan direction is used for 
several reasons. First, patterns produced are symmetric and are very easy to analyze. 
A shift off of this scan direction produces very asymmetric patterns; these patterns are 
difficult to analyze. Second, this scan direction produces patterns that can be intuitively 
compared to a uniform rectangular array steered to broadside. Finally, it is easier to 
control the array output for method B penalty functions with the symmetric patterns. 
Since the two methods investigated have different behavior and characteristics, they will 


be addressed separately. 


3.2.1 Method A 


The first parameter to be resolved for both methods is the selection of Amin or more 
specifically, the selection of the dynamic range of PF. From section 2.4.1, the larger 
Amin can be correlated with larger values of D The more the D term dominates in 
equation 2.32, the closer the penalty function resembles a sensor noise only case. This 
leads to the conventional beam pattern. Intuitively, the most interesting cases are for 
small Amin or a large dynamic range. Experience gained from working with method A 
penalty functions shows that this method has very predictable performance with respect 
to dynamic range. For the penalty function parameters investigated with method A, the 


dynamic range which leads to a sensitivity ratio of approximately 2 is on the order of 104 
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minimum eigenvalue (Amin) | dynamic range 





Table 3.2: Minimum Eigenvalues and Associated Dynamic Ranges Tested for Method A 
Penalty Functions 


with a corresponding Amin on the order of 107%, For purposes of comparison, Amin for 
method A is restricted to the values listed in table 3.2. 

With a dynamic range established, 15 cases are used to explore the effects of transition 
and center region widths on method A penalty functions. These 15 cases are summarized 
in table 3.3. Cases 1A through 5A examine the effects of center region width in azimuth. 
The best results in terms of DI and first sidelobe level occur with case 2A. For this case the 
first sidelobe level in azimuth is approximately -24.9 dB and in elevation is approximately 
-22.8 dB with DI = 37.9 and n = 2.09. For all cases 1A through 5A, changing Ad. 
has no effect on the first sidelobe level in elevation. Cases 6A through 11A examine the 
effects of center region width in elevation. Widths beyond 52.5° cause the first sidelobe 
to disappear. This results in wider main beam widths; therefore, cases beyond 52.5° are 
not considered. The best results occur with case 11A. The resulting first sidelobe level 
in azimuth is approximately -32.9 dB and in elevation is approximately -36.2 dB with 
DI = 37.4 and n = 1.62. The eigenvalue performance characteristics for case 11A in 
terms of 7 and DI are shown in figure 3-4. As opposed to cases 1A through 5A, changing - 
A0, significantly affects first sidelobe levels in azimuth where up to 8 dB differences are 
observed. This implies a partial coupling between azimuth and elevation first sidelobe 
levels. Using case 11A as the comparison point, cases 12A and 13A examine the effects 
of transition region width in azimuth while cases 14A and 15A examine the effects of 
transition region width in elevation. Transition region width appears to serve as a fine 
tuner in defining an effective method A penalty function width. The best results occur 
for case 11A; therefore, the array output, beam pattern, square magnitude and corrected 


phase for the elements weights are shown in figures 3-5 through 3-8. Figure 3-5 is the array 
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Table 3.3: Method A Case Study Summary 


processor’s representation of the penalty function. Figure 3-6 has significantly reduced 
sidelobes compared to the conventional beam pattern shown in figure 3-2. The weights, 
as defined in figures 3-7 and 3-8, do not possess a structure that is readily identified with 
a conventional weighting. 

Based on the method A performance results, the penalty function design process for 


method A is summarized by the following rules: 
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Figure 3-4: Eigenvalue Performance Results for Case 11A 
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Figure 3-5 
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Figure 3-6 
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Figure 3-7 
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e Selection of center region width is the most critical parameter choice and should be 


on the order of the width of the main beam of the conventional beam pattern. 


e Center region widths in directions with a large aperture should be selected first since 


these widths appear to be uncoupled with directions with short apertures. 


e The transition region widths should be selected to fine tune the effective penalty 
function width to yield a smooth transition over a region less than the width of the 


center region. 


e Select a dynamic range that prevents superdirective performance. Large dynamic 
ranges (but still not superdirective) usually result in lower first sidelobe levels, wider 
main beams, and smaller DI, so some trade off may be desired in the selection of 


dynamic range. 


Using these guidelines, a significant reduction in sidelobe levels can be achieved using the 


method A penalty function design process. 


3.2.2 Method B 


As was the case for method A, the dynamic range for method B penalty functions must 
first be established. This issue is not as easily addressed for method B as it is for method 
A. Although the most interesting cases are still for small Amn, method B penalty functions 
exhibit unpredictable behavior in such a way that dynamic range must be considered as 
a separate issue for each case. Therefore to cover the majority of possibilities, the values 
of Amin listed in table 3.4 are considered. These values of Amin correspond to dynamic - 
ranges varying from 10 to 10!. 

The effects of transition and center region widths for method B penalty functions 
are explored in 13 cases which are summarized in table 3.5. Cases 4B and 5B have 
performance characteristics in terms of n and DI versus Amin similar to method A; all 
are superdirective at large values of dynamic range. The remaining cases have different 
performance characteristics. These cases can not be driven superdirective even for very 
large dynamic ranges. In fact, these cases exhibit a stable performance region where no 


gain in sidelobe reduction is achieved with increasing dynamic range. The performance 
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Table 3.4: Minimum Eigenvalues and Associated Dynamic Ranges Tested for Method B 
Penalty Functions 


results for case 9B (this case displays the best results of all 13 cases) are shown in figure 
3-9. Observe the stable operating region for Amin < 107°. 

Cases 1B through 5B examine the effects of A$,. For cases 4B and 5B, A4, is too 
wide causing the first sidelobes to be emphasized in addition to the main lobe. This is 
the cause of the superdirective performance results. Case 3B results in the best sidelobe 
control with first sidelobe level in azimuth equal to -49.7 dB, first sidelobe level in elevation 
equal to -17.3 dB, Amin = 107%, n = 1.6, and DI = 34.6 dB. No effects are observed on 
sidelobe levels in elevation for cases 1B through 5B. 

Cases 6B through 8B examine the effects of A0¿. All cases produce first sidelobe levels 
within 1 dB of each other, the best being case 7B with a first sidelobe level in elevation of 
-18 dB. Sidelobe levels in elevation are generally insensitive to method B penalty function 


parameter adjustment. 
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Figure 3-9: Eigenvalue Performance Results for Case 9B 
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Table 3.5: Method B Case Study Summary 


Case 9B represents the best results of the center region width study and is used as the 
comparison case for the transition region width study. Cases 10B through 13B investigate 
the transition region width issue. Transition region effects for method B are similar to 
those effects for method A; transition region width serves as fine tuning for the effective 
penalty function width. 

Case 9B represents the best results for all 15 cases. The array output, beam pattern, 
square magnitude and corrected phase of the element weights are shown in figures 3-10 
through 3-13. Figure 3-10 shows the array processor’s interpretation of the method B 
penalty function. The beam pattern in figure 3-11 has significantly reduced sidelobes in - 
azimuth compared to the conventional beam pattern found in figure 3-2, but has only 
minor improvements in sidelobe levels in elevation. The element weights shown in figures 
3-12 and 3-13 have a structure similar to conventional weights: high in the middle and 
tapering towards the ends. 

The method B investigation also revealed several important issues characteristic of the 
method B penalty function. One, the beamformer tends to emphasize edges; in particular, 


the edges corresponding to the transition region. If the center region is wide enough, the 
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edges are located in the vicinity of the first sidelobe of the conventional beam pattern. 
This is shown in figure 3-14, the case 5B array output. This results in the superdirective 
beam pattern found in figure 3-15. More significantly, the element weights appear to have 
a threshold where the structure changes suddenly to a significantly different shape. This 
is shown in figures 3-16 and 3-17. 

The penalty function design process for method B is summarized by the following 


rules. 


e Select a center region width on the order of the width of the main beam of the 


conventional beam pattern. 


e Selection of the center region width involves a trade off against the array aperture. 
For those directions with large apertures (azimuth), the full main beam width should 
be used. For those directions with short apertures (elevation), the center region 
width must be limited to something less than the main beam width due to the wide 


main beam and low sensitivity to parameter adjustment. 


e Select transition region widths to fine tune the effective penalty function width and 


to yield a smooth transition over a region less than the width of the center region. 


e Select the smallest possible dynamic range that is located in a stable operating region 
if a stable region exists. Otherwise, choose the largest dynamic range consistent with 


superdirectivity and DI considerations. 


Using these guidelines for method B designs, a significant reduction in levels can be 


achieved with the more notable results occurring in directions of large array aperture. 


3.3 Steered Results 


The problem encountered with steered penalty functions is the introduction of asym- 
metry into the problem. This asymmetry has very little effect on method A and has 
very significant effects on method B. Only limited results are presented here; a single case 
steered to 0° elevation and 102° azimuth is illustrated. The intent is to only provide an 


intuitive feel for steered penalty function performance. 
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Figure 3-15 
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Figure 3-17: corrected phase of element weight 
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3.3.1 Method A 


Virtually no effects due to steers are encountered. When the case 9A penalty function 
is steered to 0° elevation and 102° azimuth, no distortion of the penalty function in terms 
of the array output is observed. This is shown in figure 3-18. The really nice feature of the 
method A procedure is that steering the penalty function has essentially no effect on the 
sidelobe levels compared to the broadside case. This is shown in figure 3-19. Again, this 
beam pattern has significantly better sidelobe levels than the steered conventional beam 
pattern found in figure 3-3. The square magnitude and corrected phase of the element 
weights for the steered case 9A are shown In figures 3-20 and 3-21. Again, the aperture 


weights are not characteristic of any traditional shading. 


3.3.2 Method B 


The effects of asymmetry on the method B penalty function are manifested in the 
concept of edge detections. The asymmetry causes the edges in directions of shorter 
apertures to be more emphasized which causes a significant sidelobe to be placed in the 
location of the short aperture edge. This is observed when case 9B is steered to 0° elevation 
and 102° azimuth; the array is superdirective in this case. As the penalty function center 
region 1s made wider and as the array is steered further off broadside, the asymmetry 
effects become more pronounced. (If the center region is too wide, a threshold effect 
is again observed where the weights suddenly change character into a dumbbell shape.) 
The placement of significant sidelobes at the locations of edges is similar to that which is 
observed for the broadside case 5B. This effect is illustrated by observing the array output 
and beam patterns for case 5B in figures 3-14 and 3-15 where both edges are emphasized 
instead of just one edge in the direction of short aperture as is found in the steered case 
9B. 

Based on these observations, the method B penalty function procedure for steered 
patterns must be modified. The key is to reduce the center region width to a point where 
the beamformer cannot distinguish one edge from the other; rather than two nulls in 
the array output, there appears just one center null. If Af, is reduced to 0°, then more 


reasonable results are obtained. The array output, beam pattern and 5103266 6 
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Figure 3-18 
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Figure 3-20 
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Figure 3-21 
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and corrected phase of the element weights for the zero center region width steered to 
O° elevation and 102° azimuth are shown in figures 3-22 through 3-25. Although these 
results are not as good as the broadside case, the beam pattern still has reduced sidelobes 
compared to the steered conventional beam pattern. Better performance may be obtained 


if some small center region width is inserted. 
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Chapter 4 


Summary and Conclusions 


This investigation features a systematic procedure to design beam patterns for a 
general multidimensional irregular array, an area where little previous work has been 
accomplished. The procedure uses a “penalty function” input to an MVDP beamformer. 
The proper penalty function penalizes high sidelobes and encourages the main lobe. Two 
different methods of penalty function design are investigated. Method A achieves sidelobe 
reduction by placing nulls throughout the scanning space encompassing the entire sidelobe 
region. This method is characterized by fully defining the scanning region; this leads 
to very predictable and stable performance. Method B achieves sidelobe reduction by 
emphasizing the main lobe. This is achieved by placing anti-nulls in the main beam 
region. Method B is characterized by a sparsely defined scanning region; this leads to 
sometimes unpredictable and unstable performance. The design procedure is a four step 


process: 


1. Determine the main lobe beam width for the conventional beam pattern. The main 


lobe width is best defined in terms of the angular separation between the first nulls. 


2. Select the center region width based on the main lobe beam width and the design 


method (A or B). 


3. Fine tune the effective penalty function width by selecting an appropriate transition 


region width. 
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Figure 4-1: Plot of First Sidelobe versus Center Region Width, Method A. Conventional 
main lobe widths in terms of angular separation between first nulls are approximately 24° 


in azimuth and 50° in elevation. 
4. Stabilize the beamformer by the addition of a sensor noise term in such a way to 


prevent superdirective performance. 


Selection of the center region width is the critical step. There is an optimum center region 
width with respect to the resultant first sidelobe level. This is illustrated for both methods 
as applied to the array used throughout the investigation; these results are shown in figures 
4-1 and 4-2. The figures show that a minimum in sidelobe level occurs for each direction 
except in elevation for method A where increasing the center region width eliminates the 
first sidelobe. 

Methods A and B each have distinctive characteristics. The character of each method 


must be evaluated when considering which method to apply. 


e Method A has roughly equivalent sidelobe performance in elevation and azimuth. 


Method B performance in azimuth is significantly better than performance in ele- 


vation. 


e In terms of first sidelobe level overall performance in elevation is better with method 


A; overall performance in azimuth 1s better with method B. 


e Method A results in a narrower main lobe compared to method B which results in 


better DI performance for method A. 
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Figure 4-2: Plot of First Sidelobe versus Center Region Width, Method B. Conventional 
main lobe widths in terms of angular separation between first nulls are approximately 24? 
in azimuth and 50? in elevation. 


e Method A superdirective performance is characterized by high sidelobes in the un- 


observed k space. Method B superdirective performance is characterized by high 


sidelobes in the observable k space. 
e Method A is easily steered; the method B procedure must be modified when steered. 


e For method A, azimuth performance is partially coupled to elevation performance. 


For method B, azimuth performance and elevation performance are uncoupled. 


e Method B has stable operating regions where performance is independent of penalty 


function dynamic range. 


e The method B procedure has a threshold effect where the aperture weights suddenly 
change character from shapes with one global maximum to shapes with two or more 


local maximum. 


The investigation is notably sparse of results for steered arrays; however, with the 
use of some simple cases, significant sidelobe reduction is achieved for the examined case. 
Only the effects of steers in azimuth are investigated because of the insensitivity found in 
elevation. Considering the case presented, the general requirements for steered penalty 


functions are illustrated. 
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Additional investigation into several aspects not examined are warranted. Method A 
is the design method of choice, however, it would be beneficial to obtain sidelobe reduction 
in azimuth comparable to that found in method B. One proposal to achieve these levels 
would be the incorporation of a two level background region which more heavily penalizes 
sidelobes in azimuth. 

Finally, one major issue which has been only surficially addressed concerns the robust- 
ness or sensitivity of the design procedure in terms of the uncertainty in the amplitudes 
and phases of the response. This issue is not well understood and yet is so important 
when characterizing the performance of multidimensional arrays. Further work on this 
robustness issue is required and is the next step in refining the “penalty function” design 


procedure. 
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Appendix A 


Spectral Covariance Matrix 


Inversion 


When considering an MLM beamforming approach, the problem of solving a large 
system of equations quickly arises. At the outset of this undertaking, it was hoped ex- 
ploiting the embedded Toeplitz structure of the spectral covariance matrix would lead to 
significant savings In computational time. What follows are the results of the investigation 
into Toeplitz matrix inversion techniques. 

Consider the matrix problem 


Ax=b (A.1) 


where A and b are known. The standard matrix solution method which solves any nonsin- 
gular system is Gaussian elimination.[24] The method is well documented and numerous 
algorithms are readily available which efficiently incorporate the method, particularly the | 
LINPACK routines.[12] The problem with this method is that it is generally the slowest 
of the standard methods. For an n by n real matrix, the number of multiplications for the 


Gaussian elimination method is proportional to n? 


. A Toeplitz inversion method using 
a more efficient bordering method was developed by Levinson with the number of multi- 
plications proportional to n?.[8] The goal was to incorporate these savings into the more 


complicated structure of the spectral covariance matrix. 


The multidimensional beamforming problem produces a more complicated matrix 
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scanning pattern time(sec) 


aggregate pattern 1100 
Toeplitz 4950 


single sean 265 
Toeplitz 


Table A.1: Computational Times for Matrix Solution Methods 













structure as discussed in chapter 2.2.2. Two level Toeplitz matrix solution methods suit- 
able for this complicated matrix structure have been developed.[25,13] The method used 
is from the Toeplitz Package Users Guide.[13] The number of multiplications for a real 
matrix with this method is approximately 2M?L^. For the particular array structure, 
M = 10 and L — 20. This represents a factor of 10 savings from the Gaussian elimination 
method with n — 200. Even with the overhead resulting from the permutations required 
to obtain the proper structure for input into the Toeplitz algorithm, a significant savings 
is expected. Two factors have not been accounted for up to this point. One, the matrices 
are complex. Two, the system of equations needs to be solved 1593 times. 

The two methods were tested against each other with a 1600 scanning direction pat- 
tern, where 1600 vice 1593 scanning directions were used to accommodate permutations 
of the right hand side of equation A.1, consistent with the permutation of A. The tim- 
ing results were obtained on a Digital Microvax II computer system. The run times for 
all scanning directions and for just one scanning direction are summarized in table A.1. 
Although the actual beamforming was performed on an Alliant F X40 computer system, 
the time performances for both methods were similar. 

The superior performance of the Gaussian elimination method for the full problem can 
be attributed to the LU decompositon method used. Considering the single run times, 
the full Toeplitz procedure is dependent on b and must be repeated for each new right 
hand side vector. Only a portion of the Gaussian elimination procedure is dependent on 
b. The matrix decomposition need only be accomplished once while the back substitution 


using the different right hand side vectors can be applied to the same decomposition. 


"No individual run times were examined on the Alliant, only aggregate times. The Alliant times were 
approximately 65 times faster than the Microvax. 
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In conclusion, the Toeplitz method is more efficient when the matrix problem only 
requires a few right hand side vectors; Gaussian elimination is more efficient when a large 


number of right hand sides must be evaluated. 
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Appendix B 


Dolph-Chebyshev Weights 


Aperture weights have been presented for method A and method B penalty functions. 
As previously noted, method A aperture welghts have no apparent special structure, but 
method B aperture weights have a structure which may be compared to a conventional 
welght. Since the MVDP beamformer is optimum in the minimum variance sense, a 
useful comparison might be made with Dolph-Chebyshev weights which are optimum in 
the relationship between minor lobe level and main lobe width. Using the procedures 
outlined in reference [26], the aperture weighting for a uniform rectangular array can be 
calculated with input parameters of side-lobe height in both azimuth and elevation. The 
results for -60 dB side-lobes in azimuth and -22 dB side-lobes in elevation are shown in 
figure B-1. The Dolph-Chebyshev aperture weights in figure B-1 are remarkably similar 
to the aperture weights for case 9B in figure 3-12, particularly away from the edges of 
the array. Furthermore, the side lobe levels for the case 9B beam pattern in figure 3- - 
11 can be characterized by -60 dB in azimuth and -22 dB in elevation. Based on these 
results, method B penalty function aperture weights can be favorably compared to Dolph- 


Chebyshev aperture weights. 
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